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ABSTRACT 

This  research  seeks  to  develop  non-invasive  burn  depth  evaluation  methods  from 
non-contacting  visible  and  near-infrared  spectroscopic  measurements.  In  previous  years, 
we  demonstrated  that  features  of  the  optical  reflection  spectra  of  burn  wounds  can  be 
correlated  with  the  depth  of  burn.  An  imaging  system  was  built  which  determined,  with 
accuracy  equal  to  or  better  than  that  of  a  skilled  burn  surgeon,  the  probability  that  burn 
sites  would  heal  within  three  weeks  from  date  of  injury.  Our  goal  for  the  current  project 
is  to  investigate  the  optical  reflectance  properties  of  burns,  utilizing  the  techniques  of 
multivariate  analysis,  in  order  to  improve  the  reliability  of  this  instrument. 

Excellent  progress  has  been  made  toward  achievement  of  the  goals  of  this  project. 
In  the  first  year  a  commercial  spectrophotometer  was  purchased  and  modified  to  make 
optical  reflection  measurements  in  both  the  visible  and  near-infrared  regions  (450-1800 
nm).  A  library  of  reference  spectra  was  acquired  with  this  instrument. 

During  the  second  year,  experiments  were  conducted  to  understand  the  major 
components  of  the  reflectance  spectra  of  human  skin  in  vivo.  Two  dynamic  processes 
which  influence  in  vivo  spectra  were  studied  to  improve  our  knowledge  of  burn 
physiology:  temperature  and  ischemia.  These  were  modeled  in  healthy  human  subjects. 
The  spectrum  of  pure  water,  which  is  highly  temperature  dependent  and  dominates  the 
spectra  of  biological  tissues,  was  mathematically  analyzed  for  the  purpose  of 
compensating  for  water  absorbances  in  tissues. 

Thirty  four  patients  at  Harborview  Burn  Center  were  studied,  and  the  reflectance 
spectra  of  their  burns  analyzed  with  multivariate  statistics.  An  unexpected  absorption 
band  correlating  to  burn  depth  was  observed  at  630  nm  and  was  identified  as  arising  from 
methemoglobin.  Based  on  the  foregoing  results,  a  more  reliable  algorithm  was  developed 
for  the  imaging  system.  The  new  algorithm  incorporates  wavelengths  which  measure  this 
methemoglobin  absorption.  It  has  been  tested  by  simulating  imaging  system  responses 
from  spectrophotometer  data.  The  new  model  predicts  the  healing  potential  of  1 10  burn 
sites  with  88%  accuracy,  significantly  better  than  simulated  responses  from  the  old 
imaging  system  (79%). 
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INTRODUCTION 

Burn  injuries  affect  more  than  two  million  American  victims  each  year.  Two 
hundred  thousand  of  these  will  require  hospitalization  and  10,000  will  die.1  Skin  grafting 
has  significantly  reduced  infection,  mortality,  and  disfigurement  in  patients  with  deep 
wounds.  Shallow  wounds,  however,  produce  smoother  skin  when  they  are  allowed  to 
heal  on  their  own.  As  a  general  rule,  only  burns  which  do  not  heal  in  21  days  benefit 
from  grafting2,3.  An  accurate  diagnosis  may  be  obtained  by  waiting  for  3  weeks  before 
deciding  to  operate,  but  it  has  been  shown  that  grafting  within  the  first  few  days  produces 
better  surgical  results  and  reduces  hospital  stay’2.  This  has  been  the  driving  force  behind 
the  development  of  methods  to  predict  potential  healing  times. 

For  a  wound  to  heal  well  (and  within  21  days),  it  must  re-epithelialize  from 
remnants  of  epithelium  in  the  lining  of  ducts  and  hair  follicles.  No  one  has  yet  proposed 
a  method  to  predict  healing  by  measuring  the  abundance  of  epithelium.  Instead,  most 
estimate  the  depth  of  irreparably  damaged  tissue  with  respect  to  the  total  thickness  of  the 
dermis. 

Many  hospitals  today  use  a  technique  called  clinical  assessment  which  combines 
visual  inspection  with  simple  tests,  such  as  the  wound’s  sensitivity  to  touch4.  Data 
collection,  even  for  the  largest  and  most  variegated  wounds,  requires  less  than  15 
minutes  and  does  not  cause  excessive  pain,  neither  from  the  tests  themselves  or  from 
prolonged  exposure  of  the  wound.  However,  the  accuracy  of  this  technique  is  highly 
dependent  on  the  experience  of  the  assessor.  Other  techniques,  including 
histopathological  inspection,5,6,7  topical  application  of  dyes,8  radio  isotope  studies,9 
ultrasound,10  thermography,11 12pulse  oximetry,13  laser  doppler  flowmetry,14 15 
and  multispectral  photographic  analyzers,16,17  have  proven  to  be  slower,  more  painful, 
or  less  reliable  and  are  less  commonly  used. 
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The  Bum  Depth  Indicator 

In  1980  we  built  the  Burn  Depth  Indicator  (BDI),  a  hand-held  probe  which 
illuminated  a  3  cm2  area  of  wound  with  pulsed  light  from  three  LED  chips  and  measured 
the  reflected  light.  This  technique  proved  to  be  as  fast  as  clinical  assessment,  less 
painful,  and  significantly  more  reliable.18,19  We  developed  a  similar  instrument  in  1986, 
which  provided  spatially  resolved  burn  wound  diagnosis.  This  device,  the  Imaging  Burn 
Depth  Indicator  (IBDI),  produces  a  color-coded  image  of  a  wound.20,21 22 

Selection  of  the  wavelengths  used  by  the  BDI  (green  (550  nm),  red  (640  nm),  and 
near  infrared  (880  nm))  was  based  on  the  visible  appearance  of  burn  wounds  and  the 
discovery  that  near  infrared  light  is  reflected  more  by  deep  burns  than  shallow.’16  Van 
Liew  interpreted  the  BDI  measurements  in  terms  of  blood  volume,  blood  oxygenation,  and 
eschar  thickness.  The  results  of  his  model  agreed  well  with  both  the  measured  data  and 
the  expected  physiology.*19,23 

Spectroscopic  Oximetry  in  vivo 

Measurement  of  blood  oxygenation  by  spectral  analysis  of  reflected  or  transmitted 
light  is  a  well  established  method.  Visible  (400-700  nm)  and  short-wave  near  infrared 
(700-1100  nm)  light  have  been  employed.  Hemoglobins  have  n  -  tt*,  d  -  d*  and  f  -  f* 
electronic  transitions  in  this  region.  Oxidation  and  ligand  binding  influence  the  spectrum 
substantially  (Figure  1).  A  non-invasive  oximeter,  the  Pulsed  Oximeter,24  now  common 
in  hospitals,  uses  transmission  of  660  and  920  nm  light  through  a  finger  tip.  Its  results 
correlate  well  with  arterial  blood  02  saturation  measured  by  conventional,  invasive 
methods.25  The  success  of  this  method  has  stimulated  research  on  other  non-invasive, 
spectrophotometric  oximeters.  These  new  instruments  measure  either  transmitted  or 
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reflected  light  and  have  been  applied  to  the  head,26,27  brain,28  heart,’28  and 
skin. 28,29  These  methods  are  rarely  calibrated  to  give  absolute  values,  due  to  the 
unavailability  of  reference  data,  but  they  are  quite  useful  for  clinical  trend  monitoring  and 
kinetic  studies,  where  relative  changes  are  important. 

Near  Infrared  Spectroscopy 

Near  infrared  light  (700-2500  nm)  is  absorbed  by  other  biological  components 
besides  hemoglobins.  Most  near  infrared  (N1R)  absorbances  are  due  to  overtones  and 
combination  bands  of  vibrational  transitions,  primarily  C-H,  N-H,  and  O-H.  This  region 
has  been  used  in  the  agricultural  and  food  processing  industries  to  determine,  for 
instance,  protein,  oil,  and  moisture  content.30,31  Samples  are  frequently  analyzed  in 
solid  or  powdered  form  by  diffuse  reflectance.  The  resulting  absorption  bands,  despite 
the  high  degree  of  light  scatter  in  the  samples,  are  often  linear  with  concentration. 
However  the  spectral  bands  are  broad  and  overlapping,  making  analysis  difficult. 
Calibrations  are  therefore  based  on  empirically  observed  spectral  behavior  rather  than  on 
first  principle  considerations. 

Multivariate  Statistics 

The  key  that  has  unlocked  the  potential  of  NIR  spectroscopy  is  multivariate 
statistics,  mathematical  methods  for  analyzing  data  sets  which  include  multiple 
measurements  for  each  sample.  For  spectroscopy,  each  filter  or  each  monochromator 
setting  is  used  to  obtain  different  information,  producing  a  response  vector  for  each 
sample.  Multivariate  statistics  assists  the  analyst  to  find  an  empirical  mathematical  model 
which  will  best  reproduce  the  values  of  one  or  more  constituents,  given  the  response 
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vectors  of  a  set  of  calibration  samples.  Once  the  model  is  obtained,  it  may  be  used  to 
estimate  values  for  prediction  samples,  those  whose  constituents  are  unknown. 

Constituents  are  usually  chemical  concentrations,  but  may  also  be  physical 
properties,  e.g.  crystallinity  of  nylon  yarns,  or  abstract  properties,  e.g.  the  loaf  score  of 
flour  (the  volume,  appearance,  texture,  and  spring  of  bread  made  with  that  flour)  (31). 
The  property  needs  only  to  have  some  systematic  relation  to  the  spectra.  Multivariate 
statistics  can  be  used  to  analyze  the  spectra  of  burn  injuries  to  estimate  healing  potential, 
blood  oxygenation,  and  other  properties.  Further,  the  results  of  these  analyses  and  a 
more  precise  understanding  of  the  spectra  of  burned  skin  can  be  used  to  optimize  the 
algorithm  currently  used  by  the  IBDI. 

EXPERIMENTAL 
Human  Subjects 

Forty-three  patients  were  studied  at  the  Northwest  Burn  Center,  Harborview 
Medical  Center  in  Seattle,  Washington,  over  a  13  month  period.  One  fifth  of  the  injuries 
were  selected  by  the  research  nurse.to  give  a  representative  sample  of  typical  superficial, 
shallow  partial,  deep  partial,  and  full  thickness  burns.  The  remainder  were  chosen 
because  they  were  indeterminant  or  unusual.  The  medically  fragile,  including  elderly, 
hemophiliac,  and  mentally  ill  patients,  were  excluded.  Children  were  also  excluded. 

Patients  were  examined  on  the  third  day  post-burn  as  was  done  in  previous 
studies.*18  A  few  of  the  unusual  injuries  were  studied  on  the  fourth  day,  because  the 
patients  were  transferred  from  other  hospitals  and  were  not  available  on  the  third  day. 

Two  to  five  sites  were  chosen  by  the  research  nurse,  often  with  the 
recommendation  of  the  attending  physician.  After  routine  cleaning  and  debridement 


(removal  of  loose  dead  tissue),  these  sites  were  wrapped  in  cellophane  to  prevent  drying 
and  infection.  Large  wounds  were  rebandaged,  except  for  the  selected  sites.  The  patient 
was  then  returned  to  his  own  room  and  positioned  either  in  bed  or  in  a  chair.  The 
cellophane  was  removed  prior  to  spectral  data  acquisition. 

One  spectrum  was  acquired  from  each  site,  using  an  LT  Quantum  1200 
spectrophotometer32  which  was  modified  for  this  project.  Site  locations  were  recorded 
on  polaroid  photographs.  This  procedure  generally  took  15  to  30  minutes,  including  the 
time  required  to  position  the  patient  and  to  remove  the  equipment  afterwards.  If  the 
patient  was  willing,  IBDI  images  were  also  acquired,  generally  requiring  another  10  to  20 
minutes. 

The  injuries  were  followed  by  visual  inspection  for  one  month.  The  sites  from 
which  spectra  had  been  acquired  were  classified  according  to  outcome  into  three 
categories:  Shallow  (healed  in  less  than  21  days),  Deep  (healed  in  more  than  21  days), 
and  Unknown.  Most  sites  which  were  grafted  were  classified  as  deep.  Those  which  were 
clearly  shallow  but  were  excised  due  to  their  proximity  to  a  more  severe  injury  were 
classified  as  shallow.  If  an  injury  was  not  grafted,  it  was  inspected  by  the  research  nurse 
on  or  near  the  21st  day  to  determine  if  it  had  healed.  When  that  was  not  possible,  the 
information  was  obtained  by  phoning  the  patient.  Sites  for  which  there  was  insufficient 
information  were  classified  as  unknown  and  were  excluded  from  analysis. 

Of  the  sites  whose  outcomes  were  known,  34  were  deep  and  76  were  shallow. 
The  most  common  etiology  was  flame  (52%),  followed  by  scald  (22%),  chemical  (8%), 
flash  (8%),  contact  (6%),  and  grease  (4%).  The  mean  age  of  the  patients  was  33  years 
and  the  mean  TBSA  was  6%.  There  were  3  hispanics,  2  blacks,  and  1  oriental.  Burn 
locations  included  face,  arms,  hands,  palms,  backs,  sides,  legs,  and  feet. 
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Mathematical  Methods 

The  multivariate  mathematical  methods  needed  to  interpret  data  include 
preprocessing,  calibration,  and  classification  techniques.  The  goal  is  to  find  a  useful 
relationship  between  a  sample’s  spectrum  (r)  and  the  values  of  its  constituents  (c). 
Constituents  are  usually  chemical  concentrations,  but  may  also  be  physical  properties, 
abstract  properties,  or  classifications.  The  property  needs  only  to  have  some  systematic 
relation  to  the  spectra. 

The  notation  used  will  be  as  follows:  Matrices  will  be  written  in  bold,  upper  case 
(A,  B,  C),  row  vectors  in  bold,  lower  case  (a,  b,  c),  and  scalars  in  normal  italic  {A,  B,  C, 
a,  b,  c).  Column  vectors  will  be  written  as  the  transpose  of  row  vectors  (aT,  bT,  cT).  Each 
spectrum  will  be  referred  to  as  a  response  row  vector  r,  with  length  J,  where  J  is  the 
number  of  discrete  wavelengths  contained  in  each  spectrum.  /  spectra  grouped  into  a 
matrix  will  be  R,  with  dimensions  I  by  J.  The  elements  of  R  will  be  written  as  rYl.  The 
matrix  C,  will  be  /  by  K,  where  K  is  the  number  of  spectrally  distinguishable  constituents. 
The  elements  of  C  are  cik. 

Preprocessing 

Preprocessing  transforms  the  spectra  such  that  the  spectral  features  which 
correlate  to  the  constituents  are  maximized  and  the  remainder  are  minimized.  The  first 
step  is  usually  linearization  with  respect  to  constituent  concentration.  Normally  the 
logarithm  of  instrument  response  {R)  is  used  for  linearization: 

R  =  '  /  >o  0) 

where  /0  is  the  intensity  of  the  light  incident  on  the  sample  and  /  is  the  intensity  of  the  light 
that  returns  from  the  sample  and  is  measured  by  the  instrument.  The  sample  absorbance 
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(A)  is 

A  =  -log10(ft)  (2) 

and  is  proportional  to  analyte  concentration  (c)  under  certain  conditions.  The  Beer* 
Lambert  law  is 

A  =  a  b  c  (3) 

where  b  is  the  path  le  igth,  the  distance  the  radiation  traveled  through  the  sample,  and 
a  is  the  molar  absorptivity,  a  constant  characteristic  of  the  analyte.  For  a  more  complex 
sample,  the  sample  absorbance  (A^)  is  a  linear  sum  of  the  absorbances  of  all 
components: 

Ai  *  -log10(/tj  /  /0.j)  (4) 

/|j  =  intensity  of  radiation  returned  from  the  Ith  sample  at  the  /th  wavelength 
/0  j  =  intensity  of  incident  radiation  at  the  /,h  wavelength 

Aj  =  Ajk  =  (k^  3jk  Cik) 

Ay  =  /m  sample  absorbance  at  the  /**  wavelength 

AVik  =  absorbance  of  the  kth  component  in  the  /,h  sample  at  the  /'th  wavelength 
ajk  =  absorptivity  of  the  k^  chemical  component  at  the  /th  wavelength, 

5,  *  path  length  of  light  through  the  Ith  sample 
c,k  =  concentration  of  the  k*  analyte  in  the  Ith  sample 
It  is  assumed  that  all  of  the  radiation  travels  the  same  distance  through  the  sample  so  that 
b,  is  constant  for  all  rays.  For  diffuse  reflectance,  this  is  not  strictly  true,  but  it  is  often  a 
useful  approximation. 

Other  common  transformations  include  mean  centering,  baseline  subtraction,  and 
derivatives.  The  equation  for  mean  centering  is 

('jT)m..n  centred  =  ^  -  1/J  Z^)  (6) 

where  is  the  column  vector  of  log(l  fR)  for  the  /h  wavelength  and  all  samples.  Baseline 


subtraction  is 
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(^i)ba»line  corrected  T  '  ^i,Jb  (7) 

where  r{  is  the  spectrum  for  sample  /  and  rUb  is  the  absorbance  at  the  Jb,b  wavelength. 


An  estimation  of  the  second  derivative  using  a  finite  difference  function  is 

r  =  (fh0  -  +  f  j  /  g2  (8) 

where  t •  is  the  intensity  at  the  /th  wavelength  and  g  represents  the  gap  size  of  the 
derivative. 

Calibration 

Calibration  is  used  for  continuous  constituent  values,  such  as  chemical 
concentrations  or  physical  properties.  Calibration  models  must  be  validated  before  they 
may  be  used  for  prediction.  In  calibration,  a  model  is  developed  that  relates  the  spectra 
data  matrix  (R),  composed  of  the  spectra  vectors  rj(  to  the  constituent  matrix  (C), 
composed  of  the  elements  clk,  for  a  set  of  samples,  called  the  training  set  or  the 
calibration  set  (designated  by  the  subscript  c).  A  second  set  is  used  for  validation.  The 
model  is  used  to  predict  the  constituent  values  of  the  second  set.  Poor  prediction 
indicates  that  either  irrelevant  data,  such  as  instrument  noise  was  included  in  the 
calibration  model  (overfitting)  or  that  there  is  variance  in  the  validation  set  that  was  not 
present  in  the  training  set,  such  as  a  new  component.  Once  the  model  has  been 
validated,  it  may  be  used  to  predict  the  constituents  of  unknown  samples  (designated  by 
the  subscript  p). 

Another  method  of  validation,  leave-one-out  cross  validation,  is  useful  when  there 
are  only  a  few  spectra  available.  Instead  o»  dividing  the  set  in  two,  half  for  the  calibration 
set  and  half  for  the  prediction  set,  all  samples  except  one  are  used  for  the  calibration  set 
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and  the  remaining  sample  becomes  the  prediction  set.  This  is  repeated  /  times,  leaving 
out  and  predicting  each  spectra  once.  The  result  gives  an  estimate  of  the  accuracy  of 
the  model  that  uses  all  /  samples.  Since  the  spectra  acquired  from  different  sites  on  the 
same  wound  may  be  highly  correlated,  this  approach  was  modified:  all  spectra  from  a 
single  patient  were  left  out  and  predicted  together. 

Classical  Least  Squares 

The  calibration  model  generally  consists  of  a  vector  of  matrix  of  coefficients  (S)  that 
relates  the  constituent  matrix  to  the  spectral  matrix.  The  coefficients  are  typically 
estimated  by  least  squares  or  similar  methods.  The  most  straight  forward  methods 
calculate  S  assuming  the  Beer-Lambert  relationship. 

Classical  least  squares  (CLS)  is  a  direct  calibration  technique  based  on  the  Beer- 
Lambert  law,  which  in  matrix  form  is 

R  =  C  S  +  Er  (9) 

where  S  is  the  estimated  pure  component  spectra  (sjk  a  ajk  bj  and  ER  is  the  matrix  of 
spectral  errors.  The  S  matrix  is  determined  by 

S  =  RcCc+  (10) 

where  Rc  and  Cc  are  the  spectral  and  concentration  matrices  of  the  training  set.  C+  is 
the  Moor-Penrose  pseudoinverse:33 

C+  =  (CT  C)-1  CT  (11) 

By  solving  these  equations,  ER  is  minimized  in  a  least  squares  sense.  The  estimated 
concentrations  (cpT)|  of  all  k  components  in  the  /'th  unknown  sample  with  spectrum  (rpT)j 
is 

(CpT)i  -  (r„T)i  S*  (12) 

It  is  often  useful  to  include  constraints,  such  as  requiring  all  elements  in  S  and  (cpT)j  to 
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be  greater  than  or  equal  to  zero. 

Multiple  Linear  Regression 

When  the  sample  spectra  cannot  be  approximated  as  a  simple  sum  of  pure 
component  spectra,  either  because  some  of  the  components  are  unknown  or  because 
Beer’s  law  is  inapplicable,  indirect  calibration  in  used.  The  basic  model  is 

C  =  R  P  +  Ec  (13) 

where  Ec  is  the  error  in  the  concentration  matrix  and  P  is  the  /  by  J  matrix  of  coefficients, 
similar  to  S.  C  may  be  any  constituents,  not  just  chemical  concentrations.  The  equation 
is  solved  to  minimize  Ec,  and  prediction  is 

(CpT).  -  (fpT)i  P  (14) 

Since  estimating  P  involves  calculating  the  pseudoinverse  of  Rc,  J  must  be  small.  If  the 
spectra  originally  contain  a  large  number  of  wavelengths,  there  are  two  methods  of 
simplifying  the  spectral  data:  selecting  a  few  relevant  wavelengths  and  factoring  the 
whole  spectrum. 

Multiple  linear  regression  (MLR),  also  known  as  inverse  least  squares34  and  the 
P-matrix  approach35,  is  an  inverse  least  squares  model  that  uses  the  spectral  data  from 
a  small,  selected  set  of  wavelengths.  The  model  is. 

V  =  n  p,1  +  e„  (15) 

where  p„T  is  the  vector  of  regression  coefficients  for  the  /rlh  constituent  and  ec  is  the  error 
in  the  concentration  matrix. 

The  most  common  method  for  choosing  the  wavelengths  to  use  in  the  MLR  model 
is  to  calculate  the  correlation  coefficients  between  all  wavelengths  and  the  constituent  of 
interest.  The  wavelength  with  the  greatest  correlation  in  chosen  as  the  first  wavelength 
for  the  model.  The  variance  associated  with  the  correlation  is  removed  from  the  spectral 
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and  constituent  data  and  a  second  wavelength  is  chosen,  and  so  on.  Other  methods 
also  exist. 

Principal  Component  Regression 

Principal  component  regression  (PCR)36,37  uses  the  same  fundamental  definition 
as  MLR  but  decomposes  the  full  spectral  matrix  R  into  two  factor  matrices  T  and  P, 

R  =  T  P  +  F_r  (16) 

which  can  be  calculated  by  singular  value  decomposition: 

R  =  V  S  UT  +  Er  (17) 

where  V  is  /  by  A,  S  is  A  by  A,  and  UT  is  A  by  J,  and  A  is  the  number  of  factors  (principal 
components).  T  and  P  are  VS  and  UT,  respectively.  The  matrix  T  contains  "scores," 
linear  combinations  of  concentration  data,  and  P  contains  "loadings,"  linear  combinations 
of  spectral  data.  Scores  and  loadings  are  difficult  to  interpret  in  terms  of  the  individual 
components.  The  concentrations  can  be  estimated  from  the  scores, 

C  =  T  q  +  ec  (18) 

c„  =  R„  PT  V  cc  (19) 

where  q  is  the  vector  of  coefficients  that  relates  c  to  T  and  is  estimated  as  fc+  cc.  If  cp 
is  not  calculated,  this  method  is  known  as  principal  component  analysis  (PCA)  and  does 
not  require  any  constituent  values.  This  method  is  useful  as  a  classification  technique, 
particularly  when  constituent  values  are  unavailable. 

Classification 

When  constituent  data  is  discrete,  such  as  the  place  of  origin  of  the  samples  or 
"acceptable"  and  "unacceptable",  classification  algorithms  are  used.  The  most  common 
is  K  Nearest  Neighbors  (KNN).  The  samples  are  treated  as  points  in  a  multidimensional 
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space  and  their  "nearness"  to  each  other  are  calculated: 

Dii  =  ir,  -  rj|  (20) 

where  Dti  is  the  distance  between  samples  /  and  j,  and  r,  and  ^  are  the  spectral  or  score 
vectors  for  the  two  samples.  It  is  assumed  that  samples  within  a  class  will  be  closer  to 
each  other  than  to  samples  in  another  class. 

Like  calibration,  KNN  classification  requires  a  training  set  of  samples  whose 
classifications  are  known.  However,  there  is  no  matrix  of  regression  coefficients,  an 
unknown  sample  is  simply  compared  to  all  samples  in  the  training  set,  one  at  a  time.  If 
K  =  1 ,  the  unknown  sample  is  assigned  to  the  class  of  the  training  set  sample  that  it  is 
closest  to.  Otherwise  the  classifications  of  the  K  nearest  training  samples  are  considered. 

The  data  used  for  KNN  is  often  scaled  (weighted)  so  that  the  variables  which  have 
the  largest  correlation  to  the  classification  are  bigger  than  those  which  are  less  relevant. 
One  common  method  of  determining  which  variables  are  relevant  is  the  variance 
weight,38  which  is  the  ratio  of  the  intercategory  variances  to  the  sum  of  the  intracategory 
variances.  The  equation  is 

W,„  =  2  *  [(I/Njzx,2  +  (1/Nb)Dt„2  -  (2 /NAJsysg  / 

[(1/N,)S(V*.,J2  +  (1/Nb)S(VXbJr,)2]  (21) 

where  Na  and  A/b  are  the  number  of  samples  in  classes  a  and  b,  x  are  the  values  of  the 
variable  in  question,  and  xa  m  is  the  mean  of  x  for  samples  in  class  a. 

RESULTS 
Case  Study 

A  case  study  is  included  here  to  illustrate  typical  results  for  shallow  and  full 
thickness  burns.  Patient  #14  had  a  small  scald  burn  on  her  foot.  Hot  water  had  run 
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down  into  her  high-top  boot,  producing  a  shallow  partial  thickness  burn  at  the  top  and 
a  deep  burn  where  the  water  pooled  at  the  bottom.  In  the  spectrum  of  the  shallow  partial 
thickness  burn  (Figure  2a)  the  oxy-  and  deoxy-  hemoglobin  bands  (400-600  nm)  are  very 
intense,  because  the  epidermis  was  absent  and  the  wound  was  inflamed.  In  the  near 
infrared  region  (700-1800  nm)  there  are  three  absorbance  bands  that  arise  from  the 
overtones  of  O-H  stretching  modes  of  water  (960, 1150,  and  1450  nm).  The  discontinuity 
in  the  spectrum  at  900  nm  is  an  instrumental  artifact,  caused  by  a  filter  change. 

The  spectrum  of  the  deep  burn  (Figure  2b)  was  markedly  different.  The  baseline 
was  considerably  higher  and  the  water  bands  in  the  near  infrared  were  less  intense.  The 
cause  of  these  differences  is  unknown.  The  560  nm  hemoglobin  bands  are  also  less 
intense,  as  there  is  less  blood  in  the  injured  skin.  All  of  the  above  had  been  expected. 
At  630  nm,  however,  the  is  an  unanticipated  band.  The  intensity  of  this  band,  evident  as 
a  shoulder  on  the  oxy-/deoxy-  hemoglobin  peak,  suggested  that  it  was  due  to  another 
form  of  hemoglobin.  The  band  was  attributed  to  the  acid  form  of  methemoglobin  which 
has  a  band  at  630  nm  (Figure  la).  This  compound  is  relatively  well  known  and  is  present 
in  normal  subjects  at  the  0.5-1 .0%  level.  Other  hemoglobin  derivatives  may  also  be 
present  in  smaller  quantities. 

Principal  Component  Models 
Raw  Spectra 

The  spectra  from  all  of  the  burn  wounds  are  shown  in  Figure  3.  Twelve  spectra, 
acquired  from  discolored  (black,  brown,  or  green;  n=10)  and  first  degree  (n=2)  burns 
contain  spectral  features  that  are  not  relevant  to  the  discrimination  of  shallow  and  deep 
burns.  They  were  not  included  in  any  of  the  calibration  models,  to  avoid  making  the 
models  more  complicated  than  necessary,  but  were  included  in  predictions  to 
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demonstrate  that  they  would  not  be  misclassified.  These  spectra  are  shown  in  Figure  3 
and  subsequent  figures  as  dotted  lines. 

The  spectral  regions  which  correlate  to  burn  depth  were  located  by  calculating  the 
variance  weights  (Equation  21)  for  the  absorbances  at  each  wavelength  (Figure  4).  There 
is  some  correlation  between  550  and  580  nm,  but  little  or  none  (variance  weight  =  2) 
anywhere  else.  This  indicates  that  the  hemoglobin  absorbances  (550-580  nm)  have  some 
correlation  to  burn  depth,  which  is  in  accord  with  one  of  the  fundamental  tenets  of  clinical 
assessment,  that  the  darker  pink  a  wound  is,  the  more  likely  it  is  to  be  shallow.  The  low 
correlation  in  other  portions  of  the  spectra  indicates  that  the  baseline  offsets  are  unrelated 
to  burn  depth  and  are  probably  caused  by  irreproducibility  in  the  source  /  sample  / 
detector  geometry. 

Since  the  major  spectral  features  do  not  correlate  strongly  to  burn  depth,  the 
spectra  were  analyzed  by  Principal  Component  Analysis  (Equations  16-19)  to  locate  minor 
features  which  might  be  more  relevant.  The  spectra  would  be  expected  to  display  at  least 
four  major  types  of  variation,  corresponding  to  baseline  offsets,  multiplicative  factors  (path 
length  and  total  blood  concentration),  oxygenation  of  the  blood,  and  methemoglobin 
concentration.  Minor  variations,  representing  minor  chemical  components  or  spectral 
non-linearities  would  also  be  expected.  The  first  four  principal  component  loadings 
(spectral  domain,  Figure  5a, b)  show  that  there  are  indeed  at  least  four  significant  types 
of  orthogonal  spectral  variance.  The  next  four  (Figure  5c, d)  model  only  a  small  amount 
of  the  total  variance.  It  is  possible,  but  unlikely  that  these  later  factors  are  significant. 

The  principal  component  scores  (concentration  domain,  Figure  6)  show  that  the 
first  two  factors,  are  dominated  by  irrelevant  effects  of  sample/instrument  geometry. 
Similar  effects  have  been  observed  in  the  diffuse  reflectance  spectra  of  other  irregularly- 
shaped  samples.39  The  third  and  fourth  principal  components,  however,  are  correlated 
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to  burn  depth.  The  variance  weights  for  the  first  eight  principal  components  are  2.3,  2.5, 
2.7,  3.2,  2.5,  2.3,  2.0,  and  2.0. 

The  distinction  between  "shallow"  and  "deep  is  a  somewhat  arbitrary  division  of  a 
continuous  variable,  the  fraction  of  dermis  destroyed.  It  might  be  expected  that  in  a 
multidimensional  space,  the  scores  for  shallow  burns  might  be  located  on  one  side  and 
those  for  deep  burns  on  the  other,  with  the  borderline  injuries  inbetween  (Figure  6b). 
Shallow  burns  would  most  likely  be  closest  to  other  shallow  burns,  etc.,  so  KNN  should 
give  reasonable  results.  The  KNN  algorithm  in  Pirouette40  successfully  classified  100 
(91%)  of  the  110  spectra,  using  4  principal  components.  The  optimal  number  of 
neighbors  to  consider  was  found  to  be  2.  Cross  validation,  leaving  out  one  patient’s 
spectra  each  time  (MATLAB  program  written  by  the  author),  successfully  classified  101 
(92%),  using  6  principal  components  and  considering  3  neighbors. 

Assuming  that  there  is  only  one  continuous  variable  which  relates  to  burn  depth, 
another  possible  method  of  classification  is  to  define  a  multidimensional  plane  between 
the  two  classes.  All  scores  on  one  side  of  the  plane  should  represent  shallow  burns,  and 
all  scores  on  the  other  side  should  represent  deep  injuries.  This  plane  may  be  found  by 
rotating  the  axer  r~  that  the  first  axis  corresponds  to  the  direction  of  the  variance  in  the 
data  which  most  ck  fly  distinguishes  between  the  two  classes.  The  discrimination  plane 
is  then  perpendicular  to  the  first  rotated  axis. 

The  algorithm  used  for  single-plane  discrimination  involved  5  steps.  The  first  step 
was  to  mean-center*39  the  scores.  Then  the  mean-centered  scores  were  range  scaled*39 
and  feature  weighted,*39  according  to  the  variance  weight.  The  scaled  data  was  then 
rotated,  two  dimensions  at  a  time,  starting  with  those  which  have  the  greatest  variance 
weights.  The  best  rotation  was  defined  as  that  which  produced  the  largest  variance 
weight  for  the  first  of  the  two  dimensions.  Once  the  best  rotation  was  found,  the 
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discrimination  plane  was  defined  perpendicular  to  the  first  axis,  located  so  that  the 
number  of  misclassified  calibration  samples  was  minimized.  The  samples  were  classified 
according  to  which  side  of  the  discrimination  plane  their  rotated  scores  were  located.  The 
MATLAB  code  for  these  programs  are  in  Appendix  A. 

The  scores,  rotated  in  the  first  four  dimensions,  are  shown  in  Figure  7.  The 
discrimination  plane  is  also  shown,  projected  as  a  line.  All  but  7  of  the  samples  are 
correctly  classified  by  this  plane.  Including  more  than  four  dimensions  does  not  improve 
classification,  since  only  the  first  four  principal  components  contain  much  relevant 
information.  Cross  validation,  using  4  principal  components  correctly  classified  102 
(93%).  The  single-plane  discrimination  is  slightly  better  than  KNN  because  only  one 
dimension  contains  relevant  information.  The  model  is  simpler  and  less  influenced  by  any 
single  calibration  sample. 

The  loadings,  rotated  in  the  same  manner  (Figure  8),  show  that  there  is  a  strong 
positive  correlation  between  the  absorbance  at  630  nm  and  the  depth  of  the  burn.  This 
is  probably  due  to  the  absorbance  of  methemoglobin  which  appears  to  be  characteristic 
of  deep  burns.  There  is  also  a  strong  negative  correlation  at  570  nm,  where 
oxyhemoglobin  absorbs.  It  is  reasonable  to  assume  that  the  more  healthy,  oxygenated 
hemoglobin  a  wound  contains,  the  more  likely  it  is  to  heal  quickly.  The  first  rotated 
loading  is  small  for  the  wavelengths  650-700  nm  indicating  that  this  region,  which  contains 
little  hemoglobin  absorbance,  is  not  useful  for  prediction,  This  further  supports  the 
previous  results  indicating  that  the  baseline  offsets  are  unrelated  to  burn  depth. 
Subsequent  rotated  loadings  (not  shown)  are  not  correlated  to  burn  depth  (variance 
weights  <  2.2). 
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Baseline  subtracted  spectra 

Since  the  baseline  offsets  are  useless,  they  can  be  removed  without  losing 
information,  approximating  the  baseline  by  the  apparent  absorbance  at  700  nm  (Equation 
7,  Jb  =  700  nm).  Among  the  resulting  spectra  (Figure  9)  the  soectral  differences  between 
the  two  classes  are  more  obvious.  The  intensity  of  the  hemoglobin  absorbance  (520-580 
nm)  is,  on  the  average,  greater  for  shallow  injuries  than  for  deep,  but  this  variable  alone 
is  insufficient  for  diagnosis.  Considering  only  the  intensity  of  the  absorbance,  what  the 
eye  sees,  only  69%  of  the  injuries,  at  best,  could  be  correctly  classified.  Clinical 
assessment  considers  a  few  other  variables  and  is  potentially  more  accurate. 

Since  the  variance  due  to  the  offsets  has  been  removed,  fewer  principal 
components  should  be  needed  to  describe  the  data  set.  The  first  four  loadings  (Figure 
10)  are  similar  to  those  from  the  raw  spectra  (Figure  5).  The  subsequent  loadings, 
however,  contain  more  noise.  The  second  and  third  principal  components  are  most 
highly  correlated  (Figure  11),  rather  than  the  third  and  fourth  (Figure  6),  and  the  variance 
weights  for  the  principal  components  are  greater,  2.5,  4.6,  2.9,  and  2.2  for  the  first  four. 
Prediction  results  from  KNN  and  single-plane  discrimination  are  not  significantly  better, 
but  the  model  is  simpler. 

Normalized  spectra 

Although  the  amount  of  blood  visible  in  a  wound  is  a  good  indication  of  its  viability, 
it  is  an  unreliable  indicator.  Edema  (inflammation)  and  hematoma  (bruising)  both  increase 
the  amount  of  blood  without  increasing  the  healing  potential.  In  this  study,  the  diagnostic 
errors  made  by  the  physicians  were  caused  either  by  injuries  that  were  outside  the  range 
of  their  experience,  such  as  an  unusual  chemical  burn,  or  by  deep  injuries  that  appeared 
pink  instead  of  white.  To  correctly  diagnose  these  wounds  it  is  helpful  to  focus  on  the 
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second  factor  which  relates  to  burn  depth,  the  methemoglobin  content. 

To  remove  the  variance  caused  by  variable  amounts  of  blood  in  the  skin,  the 
spectra  were  normalized  to  the  maximum  of  the  average  spectrum  (570  nm  was  chosen). 
The  normalization  equation  is 


(^i)nofmalizad  (0 baseline  subtracted  /  ^i,570nm  (22) 

where  ri  570nm  is  the  apparent  absorbance  of  the  /th  spectrum  at  570  nm.  The 
methemoglobin  absorbance  at  630  nm  is  clearly  visible  in  the  normalized  spectra  (Figure 
12). 


The  removal  of  the  multiplicative  factors  is  inaccurate  since  this  normalization  is 
somewhat  simplistic,  so  the  rank  of  the  spectral  data  set  does  not  decrease.  But  since 
the  data  is  simplified,  the  correlation  to  burn  depth  increases  (Figure  13).  The  loadings 
(Figure  14)  are  similar  to  those  of  the  baseline-corrected  spectra  (Figure  10),  but  the 
scores  of  the  first  two  principal  components  (Figure  15)  are  more  highly  correlated.  The 
variance  weights  for  the  first  four  principal  components  are  3.5,  5.3,  2.8,  and  2.2. 
Prediction  results  are  nearly  ident'cal  to  those  from  the  baseline-corrected  data. 


Hemoglobin  Spectra  Model 

Single  Layer  Model 

The  above  results  suggest  that  there  is  a  single  spectral  feature  or  single 
combination  of  interrelated  features  which  is  indicative  of  the  depth  of  a  burn  wound. 
Further,  in  the  normalized  spectra,  the  630  nm  absorbance,  attributed  to  the  acid  form  of 
methemoglobin,  correlated  to  burn  depth.  This  suggests  that  the  concentration  of  acidic 
methemoglobin  in  the  skin,  relative  to  the  total  blood  volume,  is  proportional  to  the 
amount  of  tissue  damage  and,  hence,  to  the  length  of  time  required  for  healing. 
Unfortunately,  there  is  no  way  to  test  this  hypothesis  directly.  Even  if  the  patients 


21 


consented  to  having  skin  samples  removed  for  analysis,  conventional  methods  of 
hemoglobin  speciation  are  not  applicable  to  solid  tissue  samples. 

One  alternative  is  to  attempt  to  reproduce  the  observed  reflectance  spectra, 
starting  from  the  spectra  of  pure  hemoglobins.  Assuming  that  (a)  the  particles  in  the  skin 
scatter  light  isotropically,  (b)  the  hemoglobin  is  uniformly  distributed,  (c)  there  is  no 
specular  reflectance  or  stray  light,  and  (d)  the  only  absorbing  compounds  in  the  skin  are 
hemoglobins,  then  the  reflectance  spectra  of  the  skin  could  be  modeled  as  a  linear  sum 
of  pure  hemoglobin  spectra,  using  Classical  Least  Squares  (Equations  9-1 2).  None  of  the 
assumptions  are  actually  true  and  the  hemoglobins  used  to  acquire  the  basis  spectra 
were  not  100%  pure,  so  the  model  is  not  expected  to  fit  the  data  precisely.  The  spectra 
of  the  hemoglobins  (the  matrix  S)  are  shown  in  Figure  16,  including  oxyhemoglobin, 
deoxyhemoglobin,  acid  methemogtobin,  base  methemoglobin,  and  a  baseline.  The 
spectra  were  fit  to  a  reflectance  spectrum  of  a  burn  wound,  and  a  new  spectrum  was 
reconstructed  from  the  estimated  hemoglobin  concentrations.  This  was  repeated  for  all 
burn  spectra. 

Two  examples  are  shown  in  Figure  17.  The  reconstructed  spectra  are  not  exactly 
the  same  as  the  measured  spectra,  since  the  assumptions  are  incorrect.  This  simplistic 
model  does,  however,  give  an  estimate  of  methemoglobin  in  a  burn,  information  which 
has  previously  been  unavailable.  The  resulting  predictions  are  less  accurate  (79%)  than 
those  from  the  empirical  PCA  models,  indicating  that  the  methemoglobin  concentration 
alone  is  insufficient  to  predict  healing  or  that  the  fitting  errors  were  too  large,  probably 
both.  However,  the  success  of  this  model  in  discrimination  does  strongly  support  the 
hypothesis  that  oximetric  information  correlating  to  burn  d~pth  can  be  obtained  from  the 
spectra. 
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Two  Layer  Model 

Of  the  assumptions  used  in  the  above  model,  two  are  approximately  true,  one 
might  be  true,  and  one  is  definitely  false.  The  hemoglobin  is  not  distributed 
homogeneously;  in  most  cases  there  is  a  layer  of  bloodless,  dead  tissue  (eschar)  at  the 
surface.  This  layer  reflects  incident  light,  reducing  the  amount  that  interacts  with  the 
blood  beneath.  The  reflectance  from  a  non-absorbing  layer  is  dependent  on  light 
scattering  and  is  a  function  of  wavelength. 

If  the  amount  of  light  of  wavelength  /,  reflected  from  the  sample  is  lr  and  the  total 
incident  light  is  /0  j,  then  the  instrument  response  (fl0b8j)  is 

=  h  /  'o,i  (23) 

«ob..j  =  ('s,  +  'r.j)  /  ('s0.i  +  'r,j )  (24) 

where  /r]  is  the  intensity  reflected  from  the  non-absorbing  eschar  layer,  1^  is  the  intensity 
of  the  light  that  passes  through  the  first  layer  and  enters  the  second,  blood-filled  layer, 
and  /s  j  is  the  intensity  of  the  light  returning  from  the  second  layer.  The  fraction  of  the 
incident  radiation  that  would  be  diffusely  reflected  from  the  second  layer,  if  there  were  no 


specular  reflectance  and  if  the  first  layer  were  removed,  is  fl8j, 

=  's.j  /  'so, i  (25) 

and  is  related  to  the  absorbance  of  the  blood  in  the  second  layer  (4sj), 

~  log(1/ftSij)  (26) 

It  can  be  estimated  from  the  instrument  response: 

*s,j  —  ('o,j  ^obs.j  ‘  'r,j)  /  ('o,j  "  'r.j)  (27) 

*,.i  =  (T’obs.j  *  Sj)  /  (1  -  S)  (28) 

Si  =  'r.i  /  'o.j  (29) 

if  Sj  is  known. 


For  the  burn  spectra,  S  is  not  known.  Instead,  S  was  approximated  by  a  second 
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order  polynomial,  a  function  of  the  wavelength  (A.): 

S*s0/(X°)  +  s,  /(*')  +  s2  f(X2)  (30) 

R s  should  be,  approximately,  the  product  of  the  transmittances  (Tk)  of  the  hemoglobin 
species  (/c)  times  their  concentrations  (ck). 

«.  *  «n  k  t„)  on 

And  therefore 

(7ob..i  *  Sj)  /  (1  -  Sj)  *  kn  (Ck  Tk)  (32) 

There  are  3+fr  variables  which  must  be  solved  simultaneously  (s0,  sv  s2,  cv  c2, ...  ck)  by 
using  a  non-linear  curve  fitting  algorithm  and  minimizing  the  spectral  fitting  error: 

error  =  (Tobt  j  -  Sj)  /  (1  -  Sj)  -  JI  (ck  Tk)  (33) 

The  Gauss-Newton  algorithm  supplied  by  Mathworks41,  gave  satisfactory  results,  under 
the  following  constraints: 

ck  £  0  (34) 

0  <  Sj  <  1  (35) 

The  reconstructed  spectrum  (Aobt)  is  then 

=  T„  (1-S,)  +  S,  (36) 


The  transmission  spectra  of  hemoglobins  and  the  functions  /(A.0),  /(A.1),  and  /(A.2) 
are  shown  in  Figure  18.  The  reconstructed  spectra  for  the  two  examples  are  shown  in 
Figures  19  and  20.  In  general,  the  fit  was  not  improved  if  S  was  assumed  to  be  0th  order 
with  respect  to  wavelength  (s,  =  0  and  s2  =  0).  When  S  was  1st  order,  there  was 
significant  improvement;  nearly  all  of  the  systematic  distortion  was  accounted  for.  Adding 
the  2nd  order  improved  some  of  the  more  unusual  spectra,  but  generally  had  little  effect 
on  file  "typical"  spectra.  The  1st  order  model  indicated  that  deep  burns  usually  have  more 
than  9%  add  methemoglobin  (Figure  21)  and  correctly  classified  90%  of  the  injuries. 


24 


Skin  Temperature  Model 

In  the  second  year,  spectra  were  acquired  from  the  forearms  of  healthy  volunteers. 
Their  skin  temperature  was  varied  from  9  to  35  °C  and  measured  with  a  thermistor  at  the 
surface.  The  spectra,  after  second  derivative  transformation,  were  analyzed  by  Multiple 
Linear  Regression  (Equations  13-15)  and  were  found  to  correlate  to  the  skin  temperature 
(correlation  coefficient  R  =  0.982),  with  a  standard  error  of  prediction  of  1.4  °C.  This 
experiment  is  described  in  the  previous  report.42  The  equation  was 

(  =  22.7  +  (1.21x10s)  rmsm  -  (1.44x10s)  fnmm  -  (3.85x10s)  (38) 

where  r"1015nm  is  the  second  derivative  (Equation  8,  gap  =  23)  at  1015  nm  and  t  is  the 
temperature  in  °C. 

The  model  for  estimating  skin  temperature  in  normal  skin  was  used  to  estimate  the 
temperature  of  the  burn  wounds.  The  results  suggest  that  deeper  burns  are  cooler  than 
shallow  burns  (Figure  23).  The  average  estimated  temperatures  were  32  ±  7  and  27  ± 
7  for  shallow  and  deep  burns  (P<0.01).  Thermography,  using  emitted  infrared  light,  has 
also  shown  that  deep  burns  are  cooler  than  shallow, but  that  such  measurements 
are  extremely  sensitive  to  cooling  caused  by  evaporation  of  surface  water.  NIR 
reflectance  should  be  less  sensitive  to  surface  effects,  since  it  measures  that  temperature 
within  the  skin  and  below  it.  The  model  used  here  is  not  optimal;  it  is  derived  from 
healthy  skin.  A  more  accurate  model  could  be  obtained  empirically,  if  the  temperatures 
of  wounds  could  be  measured  by  a  contacting  or  invasive  probe,  or  theoretically,  if  the 
spectral  non-linearities  were  better  understood. 

Although  several  features  in  the  NIR  spectra,  including  water  temperature  correlated  to 
burn  depth,  none  were  as  useful  as  the  hemoglobin  bands,  nor  did  they  improve 
prediction.  It  is  likely  that  there  is  still  much  useful  information  in  these  spectra,  but  in  the 
absence  of  physiological  reference  data,  this  information  could  only  be  extracted  by 
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rigorous  application  of  diffuse  reflectance  theory. 

BDI  and  IBDI  Simulations 

Now  it  is  possible  to  explain  the  predictive  ability  of  the  old  Burn  Depth  Indicators 
(imaging  and  non-imaging).  Both  instruments  used  baseline-subtracted,  single¬ 
wavelength  normalized  data.  The  baseline  offsets  were  approximated  by  the  instrument 
response  of  the  NIR  LED  (880  nm)  or  filter  (880-1100  nm)  and  the  remaining  data  was 
normalized  to  the  response  of  the  Green  (560  nm).  The  methemoglobin  concentration 
was  then,  inadvertently,  approximated  by  the  normalized  response  of  the  Red  (640  nm). 
The  exact  mathematics  are  slightly  different,  because  instruments  used  the  instrument 
response  (fl)  rather  that  log(1  /R),  but  the  principle  is  the  same. 

In  order  to  compare  the  prediction  results  from  our  spectrophotometer  with 
predictions  from  the  Burn  Depth  Indicator  (BDI)  and  the  Imaging  Burn  Depth  Indicator 
(IBDI)  built  previously,  the  BDI  and  IBDI  responses  were  simulated  from  our  spectra.  The 
emission  spectra  of  the  light  emitting  diodes  (LED’s)  in  the  BDI  were  simulated  as 
Gaussian-shaped  peaks  with  maxima  at  550,  640,  and  880  nm  and  band  widths  (fwhm) 
of  70  nm.*21*22  The  spectra  of  the  IBDI  filters  (Figure  24b)  were  obtained  from  the 
literature.*20  The  responses  for  the  near  infrared  LED  and  filter  are  somewhat  uncertain, 
because  our  spectral  data  between  750  and  890  nm  is  unreliable  and  was  not  used.  The 
simulated  results  (Figures  25a  and  25b),  however,  are  reasonably  close  to  those  obtained 
in  previous  work.*19*22  For  the  few  wounds  for  which  IBDI  images  were  acquired,  the 
simulated  results  agree  with  the  images.  The  burn  sites  used  in  this  study  seem  to  have 
a  slightly  wider  range,  both  for  the  Red/NIR  and  Green/NIR  ratios,  than  those  used  in 
previous  studies,  which  is  consistent  with  our  attempt  to  select  the  widest  variety  of 
injuries.  The  accuracy  of  the  BDI  and  IBDI  simulations  for  prediction  burn  depth  was 
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about  77  and  79%,  respectively.  VanLiew  and  Moore  reported  79  and  84%  accuracies 
for  different  sets  of  indeterminate  burns  for  the  two  instruments.*18'*20 

Proposed  Model  for  New  Imaging-Burn-Depth-Indicator 

These  results  show  that  the  model  used  for  the  Burn  Depth  Indicators  was 
essentially  sound  and  from  the  perspective  of  burn  physiology,  basically  correct,  though 
its  explanation  was  erroneous.  The  choice  of  wavelengths,  however,  could  be 
substantially  improved.  First,  the  filter  originally  selected  for  baseline  subtraction  (880- 
1100  nm)  measures  a  large  absorbance  band  at  970  nm  (Figure  24)  as  well  as  the 
offsets.  This  band  arises,  in  part,  from  the  water  in  and  on  the  skin.  The  water 
absorption  is  variable  but  not  well  correlated  to  burn  depth.  The  880-1100  nm  filter  also 
measures  the  absorbances  of  hemoglobins  (Figure  lb),  but  the  path  length  in  this  region 
is  much  longer  than  in  the  visible  region.  It  is  uncertain  whether  the  blood  in  the  volume 
sampled  is  representative  of  the  blood  in  the  skin  or  primarily  below  the  skin. 
Furthermore,  the  BDI  and  IBDI  simulations  produced  comparable  results,  in  spite  of  the 
fact  that  the  wavelength  ranges  used  for  baseline  correction  were  significantly  different. 
Therefore,  little  relevant  information  is  lost  and  some  interference  is  avoided  by  moving 
the  baseline-correction  filter  to  650  or  700  nm,  where  the  absorbances  are  significantly 
lower. 

Second,  the  normalization  filter  could  be  narrowed  and  shifted  to  produce  a  more 
reliable  response.  There  is  a  pronounced,  non-linear  distortion  of  the  spectra  which 
affects  the  550-580  nm  region  where  the  hemoglobin  absorbance  is  particularly  strong. 
In  the  spectra  of  whole  blood  (Figure  la),  the  two  oxyhemoglobin  absorption  bands  (535 
and  570  nm)  have  similar  intensities.  In  the  burn  spectra,  however,  the  535  nm  band  is 
decreased  relative  to  the  570  nm  band,  and  the  degree  of  distortion  is  unrelated  to  burn 
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depth.  The  distortion  is  probably  caused  by  scattering,  and  scattering  is  highly 
wavelength  dependent.  Therefore,  the  scattering  at  630  nm  is  likely  to  be  more  similar 
to  that  at  570  nm  than  535  nm,  and  the  former  should  give  a  more  accurate  normalization 
for  the  630  nm  methemoglobin  band.  Simulations  of  fixed-filter  instruments  using  various 
possible  normalization  filters  confirm  this.  The  normalization  LED  for  the  BDI  had  a  80  nm 
band  width  and  was  at  560  nm,  sampling  both  oxyhemoglobin  bands  equally.  The  filter 
in  the  IBDI  at  540  nm  was  worse,  sampling  only  the  most  distorted  portion  of  the 
oxyhemoglobin  spectrum.  A  filter  with  a  narrower  band  width,  30  nm  or  less,  and  located 
at  570  nm  should  give  a  more  reliable  normalization. 

Third,  the  methemoglobin  filter  should  be  optimized.  The  filter  used  previously  had 
a  80  nm  band  width  and  was  centered  at  640  nm.  The  methemoglobin  absorbance  peak 
is  at  630  nm,  and  for  the  normalized  spectra,  the  wavelength  that  correlates  most  to  burn 
depth  is  625  nm  (Figure  13).  The  band  width  of  the  methemoglobin  band  is 
approximately  40  nm.  A  narrower  filter,  40  nm  fwhm  or  less,  centered  at  630  nm  should 
give  more  precise  estimates  of  methemoglobin  absorbance. 

One  possible  filter  selection  might  include  three  20-nm  band  pass  filters  at  650, 
570,  and  630  nm.  If  the  apparent  absorbances  (log  (1/R))  of  the  three  filters  are  A650, 
A570,  and  Aggo  respectively  (Figure  26),  then  the  burn  depth  (D)  is  estimated  by 

D  ~  (^630  '  ^65o)  /  (^570  "  ^65o) 

The  discrimination  point,  calculated  the  same  way  as  the  discrimination  plane  for  the 
rotated  PCA  scores,  was  found  to  be  0.104.  Injuries  for  which  D  is  greater  than  0.104 
were  deep  and  the  rest  were  shallow.  One  hundred  one  (92%)  of  the  sites  were  correctly 
classified  by  this  model,  and  97  (88%)  were  classified  during  cross  validation.  This  model 
was  significantly  more  accurate  (88%)  than  the  IBDI  simulation  (79%)  for  the  104  sites 
simulated  (P<0.06  for  paired  data).  Small  variations  in  the  band  width  or  mean 
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wavelength  of  the  filters  do  not  significantly  affect  prediction  ability,  indicating  that  an 
instrument  using  this  algorithm  would  be  relatively  tolerant  of  manufacturing  variations. 

Although  binary  predictions  (shallow  or  deep)  are  convenient,  it  is  often  useful  to 
rank  injuries  on  a  continuous  scale,  so  that  the  size  of  the  burn  and  patient  comfort  often 
may  be  weighted  against  the  estimated  healing  potential.  The  proposed  algorithm  can 
be  used  to  produces  a  continuous  scale  (0)  as  well  as  binary  predictions.  For  sites 
which  were  known  to  have  healed  in  the  third  and  fourth  weeks  post  burn,  there  is 
roughly  a  linear  correlation  between  instrument  response  and  healing  time  (Figure  27). 
The  estimated  number  of  days  required  for  healing  (H)  is 

H  =  30  D  +  15  (40) 

This  relationship  is  not  applicable  to  sites  which  heal  in  less  than  two  or  more  than  four 
weeks,  but  such  information  would  rarely  be  important  in  management  decisions. 

For  sites  which  heal  between  14  and  28  days,  the  error  in  this  model  (±  3  days) 
is  comparable  to  the  uncertainty  of  the  reference  method  (±  2  days).  Whether  or  not  a 
wound  was  healed  (entirely  reepethilialized)  is  difficult  to  determine  precisely.  Even 
among  the  Burn  Center  staff,  opinions  differed  by  a  day  or  two.  For  wounds  whose 
outcomes  were  described  by  patients  over  the  phone  or  extrapolated  from  the  nurses’ 
last  observations,  the  uncertainty  was  about  ±3  days. 

CONCLUSIONS 

With  the  LT  spectrometer  we  are  able  to  observe  absorbance  bands  from  most  of 
the  major  constituents  of  skin.  The  instrument  is  sensitive  enough  to  monitor  changes 
as  small  as  those  induced  by  only  a  few  degrees  elevation  in  surface  temperature.  The 
techniques  used  are  completely  non-invasive,  non-contacting,  and  non-destructive  and 
require  only  a  few  minutes,  making  them  suitable  for  use  on  burn  patients. 
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In  spite  of  the  non-linear  effects  of  light  scattering  and  multiple  layers,  hemoglobin 
speciation  can  be  estimated  from  visible  reflectance  measurements  of  burn  wounds. 
Further,  these  estimates  correlate  to  the  depth  of  the  injury.  Severely  elevated  levels  of 
methemoglobin  have  been  observed  in  deep  partial  thickness  and  full  thickness  burns. 
This  suggests  that  the  Imaging  Burn  Depth  Indicator  estimates  burn  depth  by  analyzing 
the  small  color  changes  caused  by  the  presence  of  methemoglobin.  The  IBDI,  initially 
designed  to  measure  deoxyhemoglobin,  can  be  significantly  improved  by  optimizing  its 
filters  to  detect  methemoglobin.  The  resulting  instrument  would  be  fast,  accurate,  non- 
invasive,  and  non-contacting  and  could  produce  a  map  of  a  wound,  color-coded  to 
indicate  healing  potential.  This  improved  instrument  could  be  used  as  a  practical  guide 
for  burn  treatment  decisions  and  may  assist  in  the  study  of  burn  physiology. 
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Appendix  A 

MATLAB  Code  for  Rotation  and  Single-Plane  Discrimination 
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function  (Xmean, Scate.Rota]  >  spinc(X,C, Options) 

%  X  »  all  data  (samples  by  variables,  mora  samples  than  variables) 
%  C  *  class  vector,  containing  I's  and  2's  only,  0's  for  ignored  data 
%OUTUNE 
%  mean  center 
%  scale 
%  rotate 

%SUBROUTINES  NEEDED 
%  std,  varwt,  fisherwt 
%0PT10NS 

%  1  scale  (Range,  Auto,  Variance,  Fisher,  none) 

%  2  rotation  stepe 

%  3  number  of  rotated  variables  to  calculate 
%functk>n  [Xmean, Scale. Rota]  *  spinc(X,C, Options) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%CHECK  SIZES 


[m,n]  ■  size(C); 

if  m--1  &n~«1,  error  ('Rotation:  Too  many  constituents.'),  end 
ifn>m,  C»C';  (m,n]*size(C);  end 
[np.nv]  -  size(X); 

if  np-  -m  &  nv- -m,  error('Rotation:  Incompatible  matrices'),  end 
If  nv>  «m,  X-X’;  [np.nv]  »size(X);  end 
dear  m  n 


if  any(C-»1  &C--2&C--0),  error  ('Rotation:  Bad  constituent  values.'),  end 
nl  -  find(C»«l); 

if  Iength(n1)<2,  error('Rotation:  Too  few  of  class  #V),  end 
n2  -  find(C-»2); 

if  Iength(n2)<2,  error('Rotation:  Too  few  of  dass  #2'),  end 


if  ~any([1  2  3  4  5]  -  »  Options(l)), 
error(’Rotation:  Not  a  valid  scaling  option.') 
end 

if  length(OpUons)<3,  Options (3)  »nv-1;  end 
if  0ptk>ns(3)>nv-1,  Options(3)»nv-1;  end 


SETUP  MATRICES  ***** 


%  SETUP  MATRICES 
PhiStepSize  *  1  /  Optfons<2); 

Rota  -  eye(nv);  %Main  rotation  matrix 

Cycle  »  0; 


SETUP 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  PREPROCESS  DATA 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

. . »*** . * .  MEAN-CENTER  ***** 

%  MEAN-CENTER 
Xmean  «  meanCQ; 

Data  »  X-ones(np,l)*Xmean; 
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%  SCALE  (Sharaf.Mman, Kowalski;  "Chemome  tries',  p.193) 

if  Options(l)  »  ■  1,  % _  Range  scaling 

for  i«1:nv,  Scale (i)  «  1  /  max(abs(Oata(;,i))):  end 

elseif  Opttons(l)  »  »2,  % _  Autoscaling 

for  i«1:nv,  Scale (i)  -  1  /  std(Data(:,i));  end 

elseif  Options! t)  *  *3,  % _  Variance  weights 

for  l**1:nv,  DataTmp(:,i)*Data(:,i)/ma*(abs(Data(:,i)));  end 

[Scale]  -  varwt(OataTmp,C); 

for  i»l:nv,  Scale (i) «  Seale  (i)/max(abs(Data(:,i)));  end 

elseif  Options(l)  «  -4,  % _  Fisher  wsights 

for  l«l:nv,  OataTmp(:.i)»Data(:,i)/ma*(abs(Oata(:,i)));  end 
[Scale]  -  flsherwt(DataTmp,C); 
for  i»1:nv,  Scale(l)-Scale(i)/max(abs(Data(:,i)));  end 
elseif  Options(l)-  »5, 

%no  scaling 
end 

for  i»  1:nv,  Oata(:,i)«Data(:,i)*Scale(i);  end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ^%%%%%%%%%%%%%%% 

%  ROTATE  (two  vectors  at  a  time) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
forvarl  »  1.0ptions(3), 
for  vmr2  »  (varl  +  1);nv, 

% . TWO-VECTOR  SUB-MATRICES  ***** 

%Choose  vectors 
0ata2  »  0ata(:,[var1  ,var2)); 

Rota2  >  Rota(;,[var  1  ,var2]); 

Axis  -  sqrt(2)  •  (-1  1  -1  1]  *  max(max((0ata2||); 
%‘,,*****‘****‘,***************“**********************  ROTATE  (loop)  ***** 

%Spin  about  [0  0] 
for  phistep-  1:Options(2), 

%+♦  +  •*•  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  ♦  +  +  +  +  +  Rotation  Angle 
%Rotate  data 

phi  »  phistep  *  PhiStepSize  •  (2*pi); 
rot  »  [cos(phi),  -sin(phi);  sin(phi),  cos(phi)]; 
data  -  0eta2  *  rot; 
rota  »  Rota2  *  rot; 

%+  +++++++++++++++++++++++++++++++++++++++  Distance  Between  Classes 
%0istance  between  classes 
[d.sign]  -  vaiwt(data(:,l),C); 
if  sign  -  -  d  »  -d;  end 
dist (phistep)  -  d; 
end  %for  phistep 

. . *********  FIND  MAXIMUM  ANGLE  ***** 

%  Maximize  distance 

+  +  +  +  +  +  ■*•  +  +  +  +  +  ■*•  +  +  +  +  +  •*•  +  +  +  +  +  +  •»■  +  +  +  ■*•  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  Angle 
maxphi  «  find(  dist  «  «  max(dist) ); 

if  length(maxphi)  >  t,  maxphi  -  maxphi(l);  end 
phi  -  maxphi  •  PhiStepSize  *  (2*pi); 
rot  ■  (coa(phl),  -sin(phi);  sin(phi),  cos(phi)]; 

%‘**‘************************»,“ . ******  NEW  SUB-MATRICES  ***** 

%New  sub-matrices 

%+-n-  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  *-  +  +  +  +  +  +  +  New  Matrices  +  +  +  +  + 
data  «  Data2  *  rot; 
rota  »  Rota2  *  rot; 

Data(:,(var1,var2])  -  data; 

Rota(:,(var1,var2])  «  rota; 

disp(['Rotatlon:  Finished  variables  ',int2str(var1),'  and  Mnt2str(var2)]) 
end  %for  var2  ■  (varl  +  1):nv 
end  %tar  varl  -  1:(nv-1) 


function  (Data New]  *  *pinp(Xmean, Scale, Rota, Xp) 
^function  [DataNew]  -  8pinp(Xmean.Scale,Rota,Xp) 

[np  ,nv]  -  «ize(Xp); 

Data  -  Xp  -  onee(np,l)*Xmean; 

for  l»1:nv,  Data(:,i)«Data(:,i)*ScaJe(i);  end 

Data  -  Data  *  Rota; 

for  i « 1  :nv,  Data(:,i)-Data(:,i)/Scale(i);  end 
DataNew  »  Data  +  ooee(np,1)*Xmean; 
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function  [disc Point)  «  disc(0,C); 

%  CMSCc.m 

%  Finds  the  most  discriminating  point  to  separate  two  classes 
%  M.R.Balkenhol,  1/12/92 
% 

%  0  *  data  in  two  columns,  column  1  is  for  discriminating 
%  C  *  class  vector,  1's  and  2's  only 
% 

%  minE  «  best  point  (Index  and  data  value) 

%  E  »  number  of  misdassifications  [index,  classl ,  class2,  both  classes] 
% 

%function  [discPoint]  *  disc(0,C); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SETUP 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
minO  »  min(D(:,1)); 
maxD  »  max(D(:,1)); 
maxi  «  300; 

Error  -  999*ones(maxi  +  1,1); 

Error2  »  999*ones(maxi  + 1,1); 

Index  -  zeros  (maxi  + 1,1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  LOCATION  OF  MOST  DISCRIMINATING  POINT 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  i«1:maxi, 

|  »  i*(maxD-minD)/maxi  +  minO; 
el  -  sum(  D(:,1)>-j  &  C=  -1 ); 
e2  -  sum(  D(:,1)<  )  &  C»  -2  ); 

ErrorO)  «  el  ♦  e2; 
end 

[minEs]  »  min(Error); 
minis  -  find  (Error- -minEs); 


l-min(minls) :  max  (minis); 
for  i-l, 

J  «  i*(maxD-minO)/maxi  +  minD; 
nl  -  flnd(D(:,1)>-j  &C--1 ); 
n2  -  «nd(  D(:,1)<  j  &C--2); 

Error2(i)  -  sum(  (  D([n1;n2],1)  -  j  ).‘2  ); 
end 


Smallest  Sum-of-Squares  Error  ***** 


[minE, mini]  -  min(Error2); 
n  «  length(minE); 

if  n  >  1,  minE  -  minE(  round(n/2), ; );  end 
discPoint  -  minl*(maxD-minD)/maxi  +  minD; 


INPUT  and 
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function  [Cast, op]  »  discp(disc,Dp,Cp); 

%  DtSCp.m 

%  Sorts  a  predictions  tot  (Ottsi  <  disc,  Qass2  >  disc) 

%  M.RBaikonhoi,  1/12/92 
% 

%  Dp  »  data  in  two  columns,  column  1  is  for  discriminating 
%  Cp  «  class  vector,  1's  and  2‘s  only  (optional) 

% 

%  disc  »  discriminating  point,  minE(1,2)  from  disc.m 
% 

%funct»n  [Cost]  -  discp(disc,Dp,Cp); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SETUP 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
(np.nv)  -  sixe(Dp); 


%%%%%%%%%%%%%%%%%«%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  SORTING  DATA  POINTS  with  raspoct  to  DISC 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  i«1:np, 
if  Dp (1,1)  <  disc, 

Cost(i)  •  1; 
olsoif  Dp(i,i)  >  disc, 

Cost(i)  -  2; 

Costfl)  -  0; 
and 
and 

Cast  -  Cast'; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  DISPLAY  (H  Cp  exists) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  nargin  >  2, 

errors  -  find(  Cast  -  ■  Cp  ); 
plot  (find  (Cast-  - 1),  Dp(Cast«  » 1,1),  'c6+',.. 
find  (Cast-  -2),  Dp  (Cast  -  >2,1), 

[1  np],  [disc  disc), 
errors,  Dp(errort,1),  'cSo') 
ap  -  length  (errors); 
and 
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Figure  lb:  Short-wave  Near-infrared  Spectra  of  Hemoglobin 
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Figure  6.  Scores  from  principal  component  analysis  of  raw  spectra,  (x)  shallow 
and  (o)  deep  injuries. 


Figure  6.  Scores  from  principal  component  analysis  of  raw  spectra,  (x)  shallow 
and  (o)  deep  injuries. 
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Figure  10.  Loadings  from  principal  component  analysis  of  baseline-subtracted 
spectra.  Eigenvectors  a)  1  and  2,  b)  3  and  4.  Eigenvectors  (-)  1  and  3  and  (...) 
2  and  4. 
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Figure  13.  Variance  weigths  of  single-wavelength  normalized  spectra.  (— ) 
variance  weights  (•  •)  mean  spectrum. 
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Figure  15.  Scores  from  principal  component  analysis  of  normalized  spectra,  (x) 
shallow  and  (o)  deep  injuries. 
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Figure  24.a  Burn  Depth  Indicator  Simulation.  (— )  a  typical  spectrum  of  a  deep 
burn,  and  (•  •  •)  the  LED  output  functions. 
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Figure  24.b  Imaging  Burn  Depth  Indicator  Simulation.  (— )  a  typical  spectrum  of 
a  deep  burn,  and  (•  •  •)  the  filter  transmittance  functions. 
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Figure  25.  Simulated  instrument  responses,  (a)  the  burn  depth  indicator,  (b)  the 
imaging  bum  depth  indicator,  (x)  are  responses  for  shallow  bums,  (o)  for  deep 
bums.  (-)  is  the  best  discriminant  line  found  by  vanLiew. 


